Prac 5 - IT de Villiers (17502292)

In [420]:
#All the necessary imports
%matplotlib inline
import pylab as pl
import numpy as np
from scipy import signal
import IPython.display

pl.rcParams['figure.figsize'] = (9,2)

def setup_plot(title, y_label, x_label):
    pl.box(False)
    pl.margins(*(pl.array(pl.margins())+0.05))
    pl.title(title)
    pl.ylabel(y_label)
    pl.xlabel(x_label)
In [421]:
def zplane(zeros, poles):
    """
    Plot the complex z-plane given zeros and poles.
    """
    zeros=np.array(zeros);
    poles=np.array(poles);
    ax = pl.gca()

    # Add unit circle and zero axes    
    unit_circle = pl.matplotlib.patches.Circle((0,0), radius=1, fill=False,
                                 color='black', ls='solid', alpha=0.6)
    ax.add_patch(unit_circle)
    pl.axvline(0, color='0.7')
    pl.axhline(0, color='0.7')

    #Rescale to a nice size
    rscale = 1.2 * np.amax(np.concatenate((abs(zeros), abs(poles), [1])))
    pl.axis('scaled')
    pl.axis([-rscale, rscale, -rscale, rscale])
    
    # Plot the poles and zeros
    polesplot = pl.plot(poles.real, poles.imag, 'x', markersize=9)
    zerosplot = pl.plot(zeros.real, zeros.imag,  'o', markersize=9, color='none', 
                             markeredgecolor=polesplot[0].get_color(), 
                             )

    #Draw overlap text
    overlap_txt = []
    def draw_overlap_text():
        for txt in overlap_txt:
            try:    txt.remove()
            except: txt.set_visible(False)
        del overlap_txt[:]

        poles_pixel_positions = ax.transData.transform(np.vstack(polesplot[0].get_data()).T)
        zeros_pixel_positions = ax.transData.transform(np.vstack(zerosplot[0].get_data()).T)

        for (zps_pixels, zps) in [(poles_pixel_positions, poles), (zeros_pixel_positions, zeros)]:
            superscript = np.ones(len(zps))
            for i in range(len(zps)):
                for j in range(i+1,len(zps)):
                    if superscript[i]!=-1:
                        if np.all(np.abs(zps_pixels[i] - zps_pixels[j]) < 0.9):
                            superscript[i]+=1;
                            superscript[j]=-1;
            for i in range(len(zps)):
                if superscript[i] > 1:
                    txt = pl.text(zps[i].real, zps[i].imag,
                            r'${}^{%d}$'%superscript[i], fontsize=20
                            )
                    overlap_txt.append(txt)
    draw_overlap_text()

    #Reset when zooming
    def on_zoom_change(axes): draw_overlap_text()
    ax.callbacks.connect('xlim_changed', on_zoom_change)
    ax.callbacks.connect('ylim_changed', on_zoom_change)
In [422]:
import os
import urllib
import scipy.io
from scipy.io import wavfile

#Download yesterday.wav from courses.ee.sun.ac.za and return it as a numpy array
def yesterday_wav():
    url = 'http://courses.ee.sun.ac.za/Stelsels_en_Seine_414/content/yesterday.wav'
    filename = os.path.split(url)[-1]
    #Download if path does not already exist
    if not os.path.isfile(filename):
        urllib.request.urlretrieve(url, filename)
    sample_frequency, signal_array = wavfile.read(filename)
    #Normalise signal and return
    signal_array = signal_array/np.max([np.max(signal_array), -np.min(signal_array)])
    return sample_frequency, signal_array
In [423]:
#Example of subplotting
pl.figure(figsize=(11,2))
pl.subplot(1, 2, 1)
setup_plot('plot 1', 'y-label', 'x-label')
pl.stem([1,2,3,4])

pl.subplot(1, 2, 2)
setup_plot('plot 2', 'y-label', 'x-label')
pl.plot([4,3,2,1]);

Question 1: All Zeros Comb filter

In [424]:
def sinewave(F,t):
    return np.sin(2*np.pi*F*t)
fs, x1 = yesterday_wav()
In [425]:
f0 = 2205
n = np.arange(len(x1))
T = 1/fs
t = n*T
x0 = sinewave(f,t)
x2 = np.abs(x0)
x = x2+x1
In [426]:
X1 = np.abs(np.fft.fft(x1))
X2 = np.abs(np.fft.fft(x2))
X = np.abs(np.fft.fft(x))

k = np.linspace(0, fs/(2*np.pi), X.size)
k2 = np.linspace(0, 1, X.size)
time_range = np.logical_and(t>=3,t<3.005) #in order to display only between 3 and 3.005
In [427]:
pl.figure(figsize=(11,2))
setup_plot('x1(t) - Yesterday - The Beatles', '|x1(t)|', 'time (s)')
pl.plot(t[time_range], x1[time_range]) #snippet between 3s and 3.005s

pl.figure(figsize=(11,2))
setup_plot('X1(w)', '|X1(w)|', 'fw')
pl.plot(k2, X1) 
Out[427]:
[<matplotlib.lines.Line2D at 0xb5fd9668>]
In [428]:
pl.figure(figsize=(11,2))
setup_plot('x2(t)', '|x2(t)|', 'time (s)')
pl.plot(t[time_range], x2[time_range]) #snippet between 3s and 3.005s

pl.figure(figsize=(11,2))
setup_plot('X2(w)', '|X2(w)|', 'fw')
pl.plot(k2, X) 
Out[428]:
[<matplotlib.lines.Line2D at 0xb57cd3c8>]
In [429]:
pl.figure(figsize=(11,2))
setup_plot('x(t)', '|x(t)|', 'time (s)')
pl.plot(t[time_range], x[time_range]) #snippet between 3s and 3.005s

pl.figure(figsize=(11,2))
setup_plot('X(w)', '|X(w)|', 'fw')
pl.plot(k2, X) 
Out[429]:
[<matplotlib.lines.Line2D at 0x73ba8f60>]
In [430]:
IPython.display.Audio(x1, rate=fs)
Out[430]:
In [431]:
IPython.display.Audio(x2, rate=fs)
Out[431]:
In [432]:
IPython.display.Audio(x, rate=fs)
Out[432]:
In [433]:
r = 1/np.sqrt(2)
a1 = [1] + [0]*10
b1 = [1] + [0]*9 + [-1]
In [434]:
print(a1)
print(b1)
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1]
In [435]:
z1, p1, k1 = signal.tf2zpk(b1, a1)
p1= []
b1, a1 = signal.zpk2tf(z1, p1, k1)
zplane(z1, p1)
In [436]:
w1, H1 = signal.freqz(b1, a1, whole=True)
H =np.abs(H1)
setup_plot('H(w)', '|H(w)|', 'w')
pl.plot(w1, H)
Out[436]:
[<matplotlib.lines.Line2D at 0xbc543c18>]
In [437]:
y = signal.lfilter(b1, a1, x) #applying the designed filter
IPython.display.Audio(y1, rate=fs)
Out[437]:

Question 2: Comb filter with resonant poles

In [438]:
a2 = [1] + [0]*9 + [-0.95]
b2 = [1] + [0]*9 + [-1]
z2, p2, k2 = signal.tf2zpk(b2, a2)
zplane(z2, p2)
print(a2)
print(b2)
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.95]
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1]
In [439]:
w2, H2 = signal.freqz(b2, a2, whole=True)
H2 = np.abs(H2)
setup_plot('H(w)', '|H(w)|', 'w')
pl.plot(w2, H2)
Out[439]:
[<matplotlib.lines.Line2D at 0xb905b668>]
In [440]:
y2 = signal.lfilter(b2, a2, x) #applying the designed filter
IPython.display.Audio(y2, rate=fs)
Out[440]:

Conclusion

Both Comb filters succees in filtering out the noise signal. The All zero comb filter does however weaken the quality of the song, as the guitar in the background becomes quite soft and a bit distorted - thus the all zero filter, succeeds in filtering the noise, but weakens the quality of the music.

The comb filter with resonant poles, also succeeds in filtering out the noise, but in the beginning of the snippet, you can still here some of the noise, until it is filtered out completely once again. Applying this filter, does not effect the quality of the music that much.

In [ ]: